function least_squares_applications()
% LEAST_SQUARES_APPLICATIONS 最小二乘应用实例
% 
% 展示最小二乘方法在各种实际问题中的应用

fprintf('=== 最小二乘应用实例 ===\n\n');

%% 应用1：多项式拟合
fprintf('应用1：多项式拟合\n');
fprintf('目标：拟合实验数据到多项式模型\n\n');

% 生成带噪声的多项式数据
x_data = linspace(0, 2, 20)';
true_coeffs = [1, -2, 3];  % 真实系数：1 - 2x + 3x^2
y_true = polyval(fliplr(true_coeffs), x_data);
noise = 0.2 * randn(size(x_data));
y_data = y_true + noise;

fprintf('真实多项式: y = 1 - 2x + 3x^2\n');
fprintf('数据点数: %d\n', length(x_data));
fprintf('噪声水平: σ = 0.2\n\n');

% 构造Vandermonde矩阵
degree = 2;
A_poly = zeros(length(x_data), degree + 1);
for i = 0:degree
    A_poly(:, i+1) = x_data.^i;
end

% 使用QR分解求解
[Q_poly, R_poly] = qr(A_poly, 0);
coeffs_fitted = R_poly \ (Q_poly' * y_data);

fprintf('拟合结果:\n');
fprintf('估计系数: [%.4f, %.4f, %.4f]\n', coeffs_fitted);
fprintf('真实系数: [%.4f, %.4f, %.4f]\n', true_coeffs);
fprintf('估计误差: [%.4f, %.4f, %.4f]\n', coeffs_fitted' - true_coeffs);

% 计算拟合质量
y_fitted = A_poly * coeffs_fitted;
residual_norm = norm(y_data - y_fitted);
r_squared = 1 - sum((y_data - y_fitted).^2) / sum((y_data - mean(y_data)).^2);

fprintf('残差范数: %.4f\n', residual_norm);
fprintf('R²决定系数: %.4f\n', r_squared);

%% 应用2：线性回归
fprintf('\n应用2：多元线性回归\n');
fprintf('目标：预测房价基于多个特征\n\n');

% 模拟房价数据
n_samples = 100;
rng(42);  % 固定随机种子

% 特征：面积、房间数、年龄
area = 50 + 150 * rand(n_samples, 1);        % 50-200平米
rooms = 1 + 4 * rand(n_samples, 1);          % 1-5房间
age = 30 * rand(n_samples, 1);               % 0-30年

% 真实关系：价格 = 2*面积 + 10*房间数 - 0.5*年龄 + 噪声
true_price = 2*area + 10*rooms - 0.5*age + 50;
price_noise = 10 * randn(n_samples, 1);
observed_price = true_price + price_noise;

% 构造设计矩阵
X_house = [ones(n_samples, 1), area, rooms, age];
feature_names = {'截距', '面积', '房间数', '年龄'};

fprintf('特征统计:\n');
fprintf('面积: %.1f ± %.1f 平米\n', mean(area), std(area));
fprintf('房间数: %.1f ± %.1f\n', mean(rooms), std(rooms));
fprintf('年龄: %.1f ± %.1f 年\n', mean(age), std(age));
fprintf('价格: %.1f ± %.1f 万元\n', mean(observed_price), std(observed_price));

% 最小二乘回归
[Q_house, R_house] = qr(X_house, 0);
beta_hat = R_house \ (Q_house' * observed_price);

fprintf('\n回归结果:\n');
true_beta = [50, 2, 10, -0.5];
for i = 1:length(beta_hat)
    fprintf('%s: %.4f (真实: %.4f)\n', feature_names{i}, beta_hat(i), true_beta(i));
end

% 预测性能
price_pred = X_house * beta_hat;
mse = mean((observed_price - price_pred).^2);
rmse = sqrt(mse);
mae = mean(abs(observed_price - price_pred));

fprintf('\n预测性能:\n');
fprintf('均方误差 (MSE): %.2f\n', mse);
fprintf('均方根误差 (RMSE): %.2f\n', rmse);
fprintf('平均绝对误差 (MAE): %.2f\n', mae);

%% 应用3：信号去噪
fprintf('\n应用3：信号去噪\n');
fprintf('目标：从噪声信号中恢复原始信号\n\n');

% 生成原始信号
t = linspace(0, 4*pi, 200)';
signal_clean = sin(t) + 0.5*sin(3*t) + 0.3*cos(5*t);
noise_signal = 0.5 * randn(size(t));
signal_noisy = signal_clean + noise_signal;

fprintf('信号长度: %d 点\n', length(t));
fprintf('信噪比: %.2f dB\n', 20*log10(std(signal_clean)/std(noise_signal)));

% 使用傅里叶基函数进行拟合
n_harmonics = 20;
A_fourier = zeros(length(t), 2*n_harmonics + 1);
A_fourier(:, 1) = ones(size(t));  % 直流分量

for k = 1:n_harmonics
    A_fourier(:, 2*k) = cos(k*t);      % 余弦分量
    A_fourier(:, 2*k+1) = sin(k*t);    % 正弦分量
end

% 最小二乘拟合
[Q_fourier, R_fourier] = qr(A_fourier, 0);
fourier_coeffs = R_fourier \ (Q_fourier' * signal_noisy);
signal_denoised = A_fourier * fourier_coeffs;

% 去噪效果评估
snr_original = 20*log10(norm(signal_clean)/norm(noise_signal));
snr_denoised = 20*log10(norm(signal_clean)/norm(signal_clean - signal_denoised));

fprintf('去噪效果:\n');
fprintf('原始信噪比: %.2f dB\n', snr_original);
fprintf('去噪后信噪比: %.2f dB\n', snr_denoised);
fprintf('改善: %.2f dB\n', snr_denoised - snr_original);

%% 应用4：图像去模糊
fprintf('\n应用4：图像去模糊（一维示例）\n');
fprintf('目标：从模糊图像恢复清晰图像\n\n');

% 生成一维"图像"信号
n_pixels = 100;
x_img = linspace(0, 10, n_pixels)';
image_true = zeros(size(x_img));
image_true(20:30) = 1;    % 矩形脉冲1
image_true(50:60) = 0.8;  % 矩形脉冲2
image_true(70:80) = 0.6;  % 矩形脉冲3

% 模糊核（高斯模糊）
blur_width = 3;
blur_kernel = exp(-((-5:5).^2)/(2*blur_width^2));
blur_kernel = blur_kernel / sum(blur_kernel);

% 构造模糊矩阵（循环卷积）
H_blur = zeros(n_pixels, n_pixels);
for i = 1:n_pixels
    for j = 1:length(blur_kernel)
        col_idx = mod(i + j - 6 - 1, n_pixels) + 1;
        H_blur(i, col_idx) = blur_kernel(j);
    end
end

% 生成模糊图像
image_blurred = H_blur * image_true;
image_noisy = image_blurred + 0.05 * randn(size(image_blurred));

fprintf('图像尺寸: %d 像素\n', n_pixels);
fprintf('模糊核宽度: %d\n', blur_width);
fprintf('噪声水平: σ = 0.05\n');

% 最小二乘去模糊（简单逆滤波）
% 注意：实际应用中需要正则化
[Q_blur, R_blur] = qr(H_blur, 0);
image_restored = R_blur \ (Q_blur' * image_noisy);

% 评估恢复质量
mse_blur = mean((image_true - image_blurred).^2);
mse_restored = mean((image_true - image_restored).^2);

fprintf('恢复效果:\n');
fprintf('模糊图像MSE: %.6f\n', mse_blur);
fprintf('恢复图像MSE: %.6f\n', mse_restored);
fprintf('改善比: %.2f\n', mse_blur / mse_restored);

%% 应用5：参数估计
fprintf('\n应用5：非线性模型的线性化参数估计\n');
fprintf('目标：估计指数衰减模型的参数\n\n');

% 指数衰减模型：y = A * exp(-λt) + B
t_exp = linspace(0, 5, 50)';
A_true = 10;
lambda_true = 0.8;
B_true = 2;

y_true_exp = A_true * exp(-lambda_true * t_exp) + B_true;
y_noisy_exp = y_true_exp + 0.5 * randn(size(y_true_exp));

fprintf('真实参数: A = %.1f, λ = %.1f, B = %.1f\n', A_true, lambda_true, B_true);

% 线性化：ln(y - B) ≈ ln(A) - λt
% 假设B已知（实际中可能需要迭代估计）
B_assumed = 2.1;  % 稍有误差的假设值
y_linearized = log(max(y_noisy_exp - B_assumed, 0.1));  % 避免负数

% 构造线性模型：ln(y-B) = β₀ + β₁*t，其中β₀ = ln(A), β₁ = -λ
X_exp = [ones(size(t_exp)), t_exp];

% 最小二乘估计
[Q_exp, R_exp] = qr(X_exp, 0);
beta_exp = R_exp \ (Q_exp' * y_linearized);

A_estimated = exp(beta_exp(1));
lambda_estimated = -beta_exp(2);

fprintf('估计参数: A = %.2f, λ = %.2f, B = %.2f (假设)\n', ...
        A_estimated, lambda_estimated, B_assumed);
fprintf('估计误差: ΔA = %.2f, Δλ = %.2f\n', ...
        A_estimated - A_true, lambda_estimated - lambda_true);

% 验证拟合质量
y_fitted_exp = A_estimated * exp(-lambda_estimated * t_exp) + B_assumed;
r_squared_exp = 1 - sum((y_noisy_exp - y_fitted_exp).^2) / ...
                    sum((y_noisy_exp - mean(y_noisy_exp)).^2);

fprintf('拟合质量 R²: %.4f\n', r_squared_exp);

%% 应用6：数据插值和外推
fprintf('\n应用6：数据插值和外推\n');
fprintf('目标：基于稀疏数据点进行插值预测\n\n');

% 稀疏观测数据
t_sparse = [0, 0.5, 1.2, 2.1, 3.0, 4.2, 5.0]';
y_sparse = [1.0, 1.8, 2.5, 2.1, 1.5, 0.8, 0.3]' + 0.1*randn(7,1);

fprintf('观测点数: %d\n', length(t_sparse));
fprintf('时间范围: [%.1f, %.1f]\n', min(t_sparse), max(t_sparse));

% 使用样条基函数
n_knots = 10;
knots = linspace(min(t_sparse), max(t_sparse), n_knots);
A_spline = zeros(length(t_sparse), n_knots);

% 简化的B样条基函数（线性插值）
for i = 1:length(t_sparse)
    for j = 1:n_knots-1
        if t_sparse(i) >= knots(j) && t_sparse(i) <= knots(j+1)
            w = (t_sparse(i) - knots(j)) / (knots(j+1) - knots(j));
            A_spline(i, j) = 1 - w;
            A_spline(i, j+1) = w;
            break;
        end
    end
end

% 处理边界情况
A_spline(A_spline == 0) = eps;  % 避免零行

% 最小二乘拟合
[Q_spline, R_spline] = qr(A_spline, 0);
spline_coeffs = R_spline \ (Q_spline' * y_sparse);

% 在密集网格上评估
t_dense = linspace(min(t_sparse)-0.5, max(t_sparse)+0.5, 100)';
A_dense = zeros(length(t_dense), n_knots);

for i = 1:length(t_dense)
    for j = 1:n_knots-1
        if t_dense(i) >= knots(j) && t_dense(i) <= knots(j+1)
            w = (t_dense(i) - knots(j)) / (knots(j+1) - knots(j));
            A_dense(i, j) = 1 - w;
            A_dense(i, j+1) = w;
            break;
        end
    end
end

y_interpolated = A_dense * spline_coeffs;

% 评估插值质量（在观测点）
y_check = A_spline * spline_coeffs;
interpolation_error = norm(y_sparse - y_check);

fprintf('插值误差: %.6f\n', interpolation_error);
fprintf('外推范围: [%.1f, %.1f]\n', min(t_dense), max(t_dense));

%% 总结
fprintf('\n=== 应用总结 ===\n');
fprintf('1. 多项式拟合：适用于平滑数据的全局建模\n');
fprintf('2. 线性回归：多变量关系建模和预测\n');
fprintf('3. 信号去噪：频域分析和滤波\n');
fprintf('4. 图像去模糊：逆问题求解（需要正则化）\n');
fprintf('5. 参数估计：非线性模型的线性化处理\n');
fprintf('6. 数据插值：稀疏数据的密集重建\n\n');

fprintf('最小二乘方法的优势：\n');
fprintf('• 数学理论完备，有最优性保证\n');
fprintf('• 计算方法成熟，数值稳定\n');
fprintf('• 适用范围广，易于扩展\n');
fprintf('• 可以处理过定系统\n\n');

end